Can we develop phase retrieval code that can work with less than Nyquist data? We'll clearly be unable to estimate high frequency abberations (high frequency zernike modes) but we should be able to estimate the lower order ones which could still be useful for characterizing commercial / core scopes with hard to change magnifications.


In [1]:
%pylab inline
import pyotf.otf


Populating the interactive namespace from numpy and matplotlib

In [ ]:
class HanserPSFSubsample(object):
    """A class defining the pupil function and its closely related methods.

    Based on the following work

    [(1) Hanser, B. M.; Gustafsson, M. G. L.; Agard, D. A.; Sedat, J. W.
    Phase-Retrieved Pupil Functions in Wide-Field Fluorescence Microscopy.
    Journal of Microscopy 2004, 216 (1), 32–48.](dx.doi.org/10.1111/j.0022-2720.2004.01393.x)
    [(2) Hanser, B. M.; Gustafsson, M. G. L.; Agard, D. A.; Sedat, J. W.
    Phase Retrieval for High-Numerical-Aperture Optical Systems.
    Optics Letters 2003, 28 (10), 801.](dx.doi.org/10.1364/OL.28.000801)
    """

    def __init__(self, wl, na, ni, res, size, **kwargs):
        """zrange : array-like
            An alternate way to specify the z range for the calculation
            must be expressed in the same units as wavelength
        """
        nyquist_res = 1 / (2 * self.na / self.wl) / 2
        self.ratio = np.ceil(res / nyquist_res).astype(int)
        self.res = res
        self.size = size
        self.nyquist = pyotf.otf.HanserPSF(wl, na, ni, res / ratio, size * ratio, **kwargs)
        
    # do a general getter/setter that redirects to wrapped class for everything
    # aside from res

    def _gen_kr(self):
        """Internal utiltiy to generate coordinate system and other internal
        parameters"""
        k = fftfreq(self.size, self.res)
        kxx, kyy = np.meshgrid(k, k)
        self._kr, self._phi = cart2pol(kyy, kxx)
        # kmag is the radius of the spherical shell of the OTF
        self._kmag = self.ni / self.wl
        # because the OTF only exists on a spherical shell we can calculate
        # a kz value for any pair of kx and ky values
        self._kz = psqrt(self._kmag**2 - self._kr**2)

    def _gen_pupil(self):
        """Generate an ideal pupil"""
        kr = self._kr
        # define the diffraction limit
        # remember we"re working with _coherent_ data _not_ intensity,
        # so drop the factor of 2
        diff_limit = self._na / self._wl
        # return a circle of intensity 1 over the ideal passband of the
        # objective make sure data is complex
        return (kr < diff_limit).astype(complex)

    def _calc_defocus(self):
        """Calculate the defocus to apply to the base pupil"""
        kz = self._kz
        return np.exp(2 * np.pi * 1j * kz *
                      self.zrange[:, np.newaxis, np.newaxis])

    def _gen_psf(self, pupil_base=None):
        """An internal utility that generates the PSF
        Kwargs
        ------
        pupil_base : ndarray
            provided so that phase retrieval algorithms can hook into this
            method.

        NOTE: that the internal state is created with fftfreq, which creates
        _unshifted_ frequences"""
        
        self._PSFa = self.nyquist.PSFa

    # Because the _attribute_changed() method sets all the internal OTFs and
    # PSFs None we can recalculate them only when needed
    @property
    def OTFa(self):
        if self._OTFa is None:
            self._OTFa = easy_fft(self.PSFa, axes=(1, 2, 3))
        return self._OTFa

    @property
    def PSFa(self):
        if self._PSFa is None:
            self._gen_psf()
        return self._PSFa

In [2]:
import scipy.ndimage as ndi

In [3]:
ndi.zoom?


Signature:
ndi.zoom(
    input,
    zoom,
    output=None,
    order=3,
    mode='constant',
    cval=0.0,
    prefilter=True,
)
Docstring:
Zoom an array.

The array is zoomed using spline interpolation of the requested order.

Parameters
----------
input : array_like
    The input array.
zoom : float or sequence
    The zoom factor along the axes. If a float, `zoom` is the same for each
    axis. If a sequence, `zoom` should contain one value for each axis.
output : array or dtype, optional
    The array in which to place the output, or the dtype of the
    returned array. By default an array of the same dtype as input
    will be created.
order : int, optional
    The order of the spline interpolation, default is 3.
    The order has to be in the range 0-5.
mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
    The `mode` parameter determines how the input array is extended
    beyond its boundaries. Default is 'constant'. Behavior for each valid
    value is as follows:

    'reflect' (`d c b a | a b c d | d c b a`)
        The input is extended by reflecting about the edge of the last
        pixel.

    'constant' (`k k k k | a b c d | k k k k`)
        The input is extended by filling all values beyond the edge with
        the same constant value, defined by the `cval` parameter.

    'nearest' (`a a a a | a b c d | d d d d`)
        The input is extended by replicating the last pixel.

    'mirror' (`d c b | a b c d | c b a`)
        The input is extended by reflecting about the center of the last
        pixel.

    'wrap' (`a b c d | a b c d | a b c d`)
        The input is extended by wrapping around to the opposite edge.
cval : scalar, optional
    Value to fill past edges of input if `mode` is 'constant'. Default
    is 0.0.
prefilter : bool, optional
    Determines if the input array is prefiltered with `spline_filter`
    before interpolation. The default is True, which will create a
    temporary `float64` array of filtered values if `order > 1`. If
    setting this to False, the output will be slightly blurred if
    `order > 1`, unless the input is prefiltered, i.e. it is the result
    of calling `spline_filter` on the original input.

Returns
-------
zoom : ndarray
    The zoomed input.

Examples
--------
>>> from scipy import ndimage, misc
>>> import matplotlib.pyplot as plt

>>> fig = plt.figure()
>>> ax1 = fig.add_subplot(121)  # left side
>>> ax2 = fig.add_subplot(122)  # right side
>>> ascent = misc.ascent()
>>> result = ndimage.zoom(ascent, 3.0)
>>> ax1.imshow(ascent)
>>> ax2.imshow(result)
>>> plt.show()

>>> print(ascent.shape)
(512, 512)

>>> print(result.shape)
(1536, 1536)
File:      ~/miniconda3/envs/sandbox/lib/python3.6/site-packages/scipy/ndimage/interpolation.py
Type:      function

In [ ]: